Small scale structures in three-dimensional magnetohydrodynamic turbulence 
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We investigate using direct numerical simulations with grids up to 1536 3 points, the rate at which 
small scales develop in a decaying three-dimensional MHD flow both for deterministic and random 
initial conditions. Parallel current and vorticity sheets form at the same spatial locations, and 
further destabilize and fold or roll-up after an initial exponential phase. At high Reynolds numbers, 
a self-similar evolution of the current and vorticity maxima is found, in which they grow as a cubic 
power of time; the flow then reaches a finite dissipation rate independent of Reynolds number. 



Magnetic fields are ubiquitous in the cosmos and play 
an important dynamical role, as in the solar wind, 
stars or the interstellar medium. Such flows have large 
Reynolds numbers and thus nonlinear mode coupling 
leads to the formation of strong intermittent structures. 
It has been observed that such extreme events in magne- 
tohydrodynamics (MHD) are more intense than for flu- 
ids; for example, wings of Probability Distribution Func- 
tions of field gradients are wider and one observes a 
stronger departure from purely self-similar linear scaling 
with the order of the anomalous exponents of structure 
functions pj. Since Reynolds numbers are high but fi- 
nite, viscosity and magnetic resistivity play a role, tearing 
mode instabilities develop and reconnection takes place. 
The question then becomes at what rate does dissipation 
occur, as the Reynolds number increases? What is the 
origin of these structures, and how fast are they formed? 

This is a long-standing problem in astrophysics, e.g. in 
the context of reconnection events in the magnctopausc, 
or of heating of solar and stellar corona. In such fluids, 
many other phenomena may have to be taken into ac- 
count, such as finite compressibility and ionization, lead- 
ing to a more complex Ohm's law with e.g. a Hall cur- 
rent or radiative or gravitational processes to name a few. 
Many aspects of the two-dimensional (2D) case are un- 
derstood, but the three-dimensional (3D) turbulent case 
remains more obscure. Pioneering works 0] show that 
the topology of the reconnecting region, more complex 
than in 2D, can lead to varied behavior. 

The criterion for discriminating between a singular and 
a regular behavior in the absence of magnetic fields fol- 
lows the seminal work by Beale, Kato and Majda (here- 
after BKM) where, for a singularity to develop in the 
Euler case, the time integral of the supremum of the vor- 
ticity must grow as (t— t t )~ a with a > 1 and the singu- 
larity time. In MHD Q , one deals with the Elsasser fields 
z ± = v±b and uj^ = w±j = Vx (v±b), with u> the vor- 
ticity, v the velocity, j the current density and b = V x A 
the induction in dimcnsionlcss Alfvcnic units, A being 
the vector potential. Intense current sheets are known to 
form at either magnetic nulls (b = 0) or when one or two 
(but not all) components of the magnetic field go to zero 
or have strong gradients. In two dimensional configura- 
tions, a vortex quadrupole is also associated with these 
structures. The occurrence of singularities in MHD has 



been examined in direct numerical simulations (DNS), 
with either regular 0, [|| or adaptive grids , and with 
different initial configurations with no clear-cut conclu- 
sions in view of the necessity for resolving a large range of 
scales (see @ and references therein for the Euler case). 
Laboratory experiments and DNS have also studied the 
ensuing acceleration of particles in the reconnection re- 
gion (see e.g. @)- 

The early development of small scales in such flows is 
exponential |lOl | (in the context of turbulent flows, see e.g. 
[ll|), because of the large-scale gradients of the velocity 
and magnetic fields, assumed given, stretching the vor- 
ticity and current. The phase beyond the linear stage, 
though, is still unclear. In 2D, numerical simulations 
with periodic boundary conditions show that the late- 
time evolution of non-dissipative MHD flows remains at 
most exponential [l2T | , a point latter confirmed theoreti- 
cally [13J by examining the structure around hy per bolic 
nulls, although finite dissipation seems to set in 

In 3D, most initial conditions develop sheets that may 
render the problem quasi two-dimensional locally; 3D 
MHD flows display a growth of small scales of an ex- 
ponential nature, although at later times a singular be- 
havior may emerge @. In this light, we address in this 
paper the early development of structures in 3D and the 
ensuing evolution in the presence of dissipation. 

The incompressible MHD equations read: 
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together with V-v = 0, Vb = 0; V is the pressure, po = 1 
is the constant density, and v and r\ are the kinematic 
viscosity and magnetic diffusivity. With v = 0, r/ = 0, 
the energy E = (y 2 + b 2 ) 1 2 and cross helicity He = 
(v • b) /2, are conserved [T3, with the magnetic helicity 



H M = (A • b) in 3D. Defining D±/Dt = d t 
can symmetrize Eqs. and obtain [Tlj : 
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TABLE I: Runs with an Orszag-Tang vortex (OT1-4), or 
with large-scale ABC flows and small-scale random noise with 
a k~ z spectrum (RND1-5); N is the linear resolution. 



Run 


N 3 


v — rj 


OT1 - OT4 


64 J - 512 3 


1 x 10~' A - 7.5 x 10" 4 


RND1 - RND4 


64 3 - 512 3 


8 x 1CT 3 - 6 x 10~ 4 


RND5 


1536 3 


2 x 1(T 4 



omitting dissipation. Note that the first term on the 
r.h.s. of (j2J) is equal to zero in 2D; the second term is 
absent in the Navier-Stokes case and may account for 
extra growth of the generalized vorticities for conducting 
fluids unless the Elsasser field gradients are parallel. 

To study the development of structures in MHD turbu- 
lence, we solve numerically Eqs. using a pseudospec- 
tral method in a three dimensional box of side 27r with 
periodic boundary conditions. All computations are de- 
aliased, using the standard 2/3 rule. With a minimum 
wavenumber of k m i n = 1 corresponding to Lq = 27r, a lin- 
ear resolution of N grid points has a maximum wavenum- 
ber k m ax = N/3. At all times we have ko/kmax < 1, 
where is the dissipation wavenumber evaluated using 
the Kolmogorov spectrum (at early times the resolution 
condition is less stringent). Two different initial condi- 
tions are used; Table [I] summarizes all the runs. 

As the system is evolved, we monitor the small scale 
development by following the dynamical evolution of the 
extrema of the generalized vorticities or of their individ- 
ual components [3 in the spirit of the BKM criterion. 

We start discussing the results for the Orszag-Tang 
vortex (OT hereafter); in two dimensions ^ has be- 
come a prototype flow for the stud y o f MHD turbulence, 
including in the compressible case In 3D, the veloc- 
ity and the magnetic field are taken to be: 

vo = [— 2 sin y, 2 sinx, 0] , 

bo = P[— 2 sin2j/ + sin z, 2 sinx + sinz, sin a; + siny]. 

The OT flow in 2D has a stagnation point in the (x, y) 
plane and an hyperbolic X-point for the magnetic field; 
a 3D perturbation is added in the z direction, resulting 
in a flow that has nulls for the magnetic field (three com- 
ponents equal to zero) of different types fljil ] correspond- 
ing to the signs of the eigenvalues of the dibj matrix p|; 
initially, the kinetic and magnetic energy Ey and Em 
are equal to 2 with f3 = 0.8, the normalized correlation 
p c = 2 (v • b) /{(v 2 + b 2 )) ~ 0.41, and H M = 0. 

Four runs were done for spatial resolutions up to 512 3 . 
The Reynolds number R e = UL/v (where U = (v 2 ) is 
the rms velocity, L = 2ir J Ev(k)k~ 1 dk/ J Ev(k)dk is 
the integral scale, and Ev(k) is the kinetic energy spec- 
trum) ranges from 570 to 5600 at the time of maximum 
dissipation of energy e = — v (w 2 ) — rj (j 2 )- FigureQJa) 
shows the temporal evolution of the maximum of the cur- 
rent max{j} (the vorticities u> and u)^ behave in a sim- 
ilar fashion). After an initial exponential phase up to 
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FIG. 1: (a) Evolution of the supremum of current for the OT 
runs in log-log. The inset shows the evolution at early times 
in lin-log units; a slope of t 3 is also indicated. The exponential 
phase ends at t ~ 0.6. (b) Total dissipation as a function of 
time for the same runs in lin-lin. i? e = 570 (solid), R e = 1040 
(dot), R e = 3040 (dash), and R e = 5600 (dash-dot). 



t ~ 0.6 and corresponding to the linear development of 
current (and vorticity) sheets through stretching by ve- 
locity gradients, a faster increase develops with, at high 
resolution, a self-similar ~ t 3 law. Note that the growth 
of max{j} during the early exponential phase seems to 
be independent of the value of R e . 

The first temporal maximum of max{j} is reached at 
slightly later times as R e increases; similarly [see Fig. 
IHb)], the total energy dissipation e shows a delay in the 
onset of the development of small scales as R e grows, rem- 
iniscent of 2D behavior [T^j , with a slower global growth 
rate after the initial exponential phase; this delay, how- 
ever, does not preclude reaching a quasi-constant maxi- 
mum of e in time as R e grows. Whereas in the 2D case, 
the constancy of e only obtains at later times when re- 
connection sets in, with a multiplicity of current sheets, 
in the 3D case more instabilities of current and vorticity 
structures are possible and the flow becomes complex as 
soon as the linear phase has ended. The dependence on 
R e of the time at which the first maximum of max{j} is 
reached is slow (~ i?° 08 ), and similarly for the time the 
maximum of e is reached. Computations at higher R e 
should be performed to confirm these results. 

The sharp transition around t = 0.6 can be interpreted 
in terms of the non-locality of interactions in MHD tur- 
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FIG. 2: (a) Evolution of the supremum of the current density 
in log-log for runs RND1 to RND5, with R e = 690 (solid), 
R e = 1300 (dot), Re = 2300 (dash), R e = 4200 (dot-dash), 
and R e = 10100 (long dash). At high R e , a power law con- 
sistent with t 3 is recovered, (b) Magnetic energy spectra at 
early times in run RND5. The lines (from below) correspond 
to t — 0.6 up to t = 1.6 with temporal increments of 0.2. 
Slopes of k~ 3 and fc~ 5//3 are indicated as a reference. 

bulence |20| with transfer of energy involving widely sep- 
arated scales. Thus, as the flow exits the linear phase, 
all scales can interact simultaneously; this may be a rea- 
son why, in a computation using OT and adaptive mesh 
refinement 0, it was found to be difficult to go beyond 
t = 0.6 since small scales were developing abruptly in 
many places in the flow through prevailing nonlincaritics. 
The energy spectra in this early phase are steep, with a 
power law ~ k~ 3 (not shown). A shallower ~ fc -170 spec- 
trum develops at later times, as found in earlier works. 

In view of similarities between the behavior observed 
on the 3D OT vortex and its 2D counter-part, it is worth 
asking whether such a development is not due to the 
high degree of symmetry of the flow. In that light, we 
now examine the temporal development of a Beltrami 
flow on which small scale random fluctuations are added. 
The initial velocity and magnetic field spectra arc taken 
~ fc-3 e -2(fc/fc ) 2 . tne s hells with k £ [1,3] have a su- 
perposition of three ABC flows |2l|, and the rest of the 
spectrum is loaded with Fourier modes with Gaussian 
random realizations for their phases chosen so that ini- 
tially, E v = E M = 0.5, p c ~ 10~ 4 and H M ~ 0.45. 
Unlike the OT runs, there are no 2D null points or ex- 
act zeros in the magnetic field. As for the OT case, four 



FIG. 3: (Color online) Regions of strong current density, and 
magnetic field lines in their vicinity for run RND5 at t = 1.6. 
The region at left has 450 2 x 250 points, and that at right has 
260 x 160 x 200. The sheets are thin and elongated (up to 1/3 
the size of the box); the magnetic field lines are parallel to 
the sheet and quasi-orthogonal to each other on each side of 
it, and they depart from the sheet transversally. Both folding 
(left) and rolling (right) occurs at this R e . Vortex sheets (not 
shown) are co-located and parallel to the current sheets. 

runs (RND1-RND4) were done with resolutions ranging 
from 64 3 to 512 3 grid points; RND5 on a grid of 1536 3 
points is run until saturation of growth of the maximum 
current. 

Both the exponential and the self-similar phases are 
noisier (see Fig. as can be expected with several 

structures competing for the small scale development of 
maxima. At low i? e , self-similar evolution seems to occur 
at a slower pace, with laws ~ t 2 , as in fact also found in 
2D at comparable resolutions. However, the two runs 
with highest resolution (RND4 and RND5) indicate a 
steeper development compatible with a t law. 

Figure|2Jb) shows the evolution of the magnetic energy 
spectrum EmOs) at early times, during the self-similar 
growth of the current density [the evolution of Ey(k) is 
similar]. Before t ~ 0.6, the largest wavenumbcrs have 
amplitudes of the order of the truncation error. For t > 
0.6, as all scales are nonlincarly excited, a self-similar 
growth sets in and the energy spectra are compatible 
with a fc -3 law. After max{j} saturates, the slope of 
Em(1c) increases slowly towards a fc -170 law. The same 
behavior is observed in the OT run, in which no fc~ 3 
power law is imposed in the initial conditions. 

The structures that develop appear to be accumula- 
tions of current sheets (similarly for the vorticity, not 
shown) , as was already found in || . Figure shows a 
zoom on two such structures, with the magnetic field lines 
indicated as well. It appears clear from such figures that 
only one component of the magnetic field reverses sign 
in most of these configurations, reminiscent of magneto- 
spheric observations. Both terms appearing in Eq. Q 
for the dynamical evolution of uj^ 1 arc substantial and 
comparable in magnitude although they may be quite 
weak elsewhere in the flow. Kclvin-Hclmoltz instabili- 
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ties with rolling up of such sheets are also present in 
the flow but only at the highest Reynolds number (run 
RND5); at lower R e the sheets are thicker, the insta- 
bility is too slow and only folding of such sheets occur. 
Magnetic field lines are parallel to the roll, in such a 
way that magnetic tension does not prevent the occur- 
rence of the instability. Note that folding of magnetic 
structures has been advocated in the context of MHD 
at large magnetic Prandtl number . Alfvenization of 
the flow (v = ±b) is rather strong in the vicinity of the 
sheets, with 0.7 < \p°\ < 1, although globally the flow 
remains uncorrelated (p c ~ 4 x 10~ 4 ); this local Alfv- 
enization gives stability to such structures since the non- 
linear terms are weakened, in much the same way vor- 
tex tubes in Navier-Stokes flows are long-lived because 
of (partial) Beltramization (v ~ ±u>). Moreover, within 
the sheet p c is positive, and it is negative outside, with a 
slight predominance of b. All this indicates that a double 
velocity-magnetic field shear plays an important role in 
the development of small scales in MHD. 

There is an elementary, analytically-soluble, one- 
dimensional model that illustrates sharply the role that 
velocity shear can play in enhancing^ current intensity, 
e.g. during early dynamo activity [23j- This consists of 
two semi-infinite slabs of rigid metal with equal conduc- 
tivities, at rest against each other at the plane y = 0, 
say. A uniform dc magnetic field bo is perpendicular to 
the interface and penetrates both slabs. At time t = 0, 
the slabs are impulsively forced to slide past each other 
in the x-direction with equal and opposite velocities (vq , 



say). The developing (quasi-magnetostatic) field, which 
acquires an x-componcnt, is a function of y and t only, 
and is governed by diffusion equations above and below 
the plane y = 0. Matching tangential components of the 
electric field immediately above and below the interface 
reduces the pair to a soluble diffusion equation with a 
fixed y-derivative at y = 0. The resulting magnetic field 
is expressible in terms of complementary error functions 
and grows without bound, as does the total Ohmic dis- 
sipation. The introduction of a time dependence in vo 
may allow also for solutions in which the maximum of 
the current grows as a power law in time. 

In conclusion, high resolution simulations of the early 
stages in the development of MHD turbulence allowed us 
to study the creation and evolution of small scale struc- 
tures in three-dimensional flows. Roll up of current and 
vortex sheets, and a self-similar growth of current and 
vorticity maxima was found, features that to the best 
of our knowledge were not observed in previous simula- 
tions at smaller Reynolds numbers. Also, a convergence 
of the maximum dissipation rate to a value independent 
of R e was found. More analysis will be carried out to 
understand how structures are formed, the relevance of 
the development of alignment between the fields and the 
creation and role of local exact solutions to the MHD 
equations (such as Alfven waves). 
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